/*
 *  composite.cpp
 *  EPPI-0.0
 *
 *  Created by Sergio Conde Martin on 6/02/13.
 *  Copyright 2007 __MyCompanyName__. All rights reserved.
 *
 */

#include "composite.h"

Composite::Composite()
:	Material(),
	va0(),
	A0(),
	nf(),
    kappa(),
    k1(),
    k2(),
	matGround()
{}

// ***

Composite::Composite(ifstream &in, const string& na, const string & t, Material* ma)
:	Material(na,t),
 	matGround(ma),
	va0(),
	k1(),
	k2(),
	kappa(1.0/3.0),
	nf(1),
	IVm(1.0),
	VIm(1.0),
	Wf(0.0),
    WfIm1(0.0),
    WfIm2(0.0),
    WfIVm(0.0),
    WfVIm(0.0),
    WfIm1Im1(0.0),
    WfIm2Im2(0.0),
    WfIm1IVm(0.0),
    WfIm2VIm(0.0),
    WfIVmIm1(0.0),
    WfVImIm2(0.0),
    WfIVmIVm(0.0),
    WfVImVIm(0.0)
{
	int i;
	string keyword;
	string line;
	vector<string>	aValues;
	stringstream stnum;
	double value;
	
	while (in >> keyword && keyword.compare("ENDMATERIAL") )
	{
        
        if(!keyword.compare("FIBERSDIR"))
		{
			getline(in,line);
			istringstream iss(line);
			while (iss >> keyword )
			{
				aValues.push_back(keyword);
			}
			nf = aValues.size() / 3;

			for(i=0; i < 3*nf; ++i)
			{
				stnum << aValues[i];
				stnum >> value;
				va0.push_back(value);
				stnum.clear();
			}
		}
        if (!keyword.compare("GHO-STRESS"))
        {
            in >> k1;
        }
        if (!keyword.compare("GHO-DLESS"))
        {
            in >> k2;
        }
        if (!keyword.compare("DISPERSION"))
        {
            in >> kappa;
        }
    }

   // double norm2a0;
    
    switch (nf)
    {
        case 1:
            
            norm2a0 = va0[0]*va0[0] + va0[1]*va0[1] +va0[2]*va0[2];
            
            A0(1,1) = va0[0] * va0[0] / norm2a0;
            A0(2,2) = va0[1] * va0[1] / norm2a0;
            A0(3,3) = va0[2] * va0[2] / norm2a0;
            A0(2,1) = va0[0] * va0[1] / norm2a0;
            A0(3,1) = va0[0] * va0[2] / norm2a0;
            A0(3,2) = va0[1] * va0[2] / norm2a0;
            
            break;
            
        case 2:
            
           // double norm2b0;
            
            norm2a0 = va0[0]*va0[0] + va0[1]*va0[1] +va0[2]*va0[2];
            norm2b0 = va0[3]*va0[3] + va0[4]*va0[4] +va0[5]*va0[5];
            
            A0(1,1) = va0[0] * va0[0] / norm2a0;
            A0(2,2) = va0[1] * va0[1] / norm2a0;
            A0(3,3) = va0[2] * va0[2] / norm2a0;
            A0(2,1) = va0[0] * va0[1] / norm2a0;
            A0(3,1) = va0[0] * va0[2] / norm2a0;
            A0(3,2) = va0[1] * va0[2] / norm2a0;
            
            B0(1,1) = va0[3] * va0[3] / norm2b0;
            B0(2,2) = va0[4] * va0[4] / norm2b0;
            B0(3,3) = va0[5] * va0[5] / norm2b0;
            B0(2,1) = va0[3] * va0[4] / norm2b0;
            B0(3,1) = va0[3] * va0[5] / norm2b0;
            B0(3,2) = va0[4] * va0[5] / norm2b0;
            
            break;
    }
    
    cout << "Number of fibers= " << nf << endl;
    cout << "Dispersion= " << kappa << endl;
    cout << "k1= " << k1 << endl;
    cout << "k2= " << k2 << endl;
    cout << "A0= " << endl << A0 << endl;

    hypertype = matGround->hypertype;

    Sf.setZero(6);
    Ctf.setZero(6,6);

    // Atributtes from the base
    Siso.setZero(6);
    Svol.setZero(6);
    Ctiso.setZero(6,6);
    Ctvol.setZero(6,6);

    density = matGround->getDensity();

    // cout << "matGround: " << matGround->type << endl;
}

// ***

void Composite::calcAnisoTerms( const VectorXd & C )
{
    // Calculation of invariants

	//matGround->calcInvariants(C);			//Comprobar si me lo puedo ahorrar

	Im = pow(matGround->getJ(), -2.0/3.0) * matGround->getInvariants()(0);
	//matGround->calcCTensor(C);
    IVm = pow(matGround->getJ(), -2.0/3.0) * doubleContraction(A0, matGround->getCTensor());
    VIm = pow(matGround->getJ(), -2.0/3.0) * doubleContraction(B0, matGround->getCTensor());

  //cout << "Invariants= " << endl << Invariants << endl;
  //cout << "Im= " << Im << endl;
  //cout << "I4= " << IVm / Jm << endl;
  //cout << "I6= " << VIm / Jm << endl;
	
	//if (Jm <= 0.0) cout << "*ERROR: Negative volume, material " << name << " has collapsed \n";
	
    // Calculation of derivatives
    
    if(IVm > 1.0)
    {
    	double exp_1 = exp( k2 * pow((1.0 - 3.0 * kappa) * IVm + kappa * Im - 1.0, 2.0));


        WfIm1 = k1 * kappa * (kappa * Im + (1.0 - 3.0 * kappa) * IVm - 1.0) * exp_1;
        WfIVm = k1 * (1.0 - 3.0 * kappa) * (kappa * Im + (1.0 - 3.0 * kappa) * IVm - 1.0) * exp_1;
        
       //cout << "WfIm1= " << WfIm1 << endl;
       //cout << "WfIVm= " << WfIVm << endl;
        
        // Calculation of second derivatives
        
        WfIm1Im1 = 2.0 * k1 * k2 * pow(kappa,2.0) * pow((1.0 - 3.0 * kappa) * IVm + kappa * Im - 1.0, 2.0) * exp_1 + k1 * pow(kappa, 2.0) * exp_1;
        
        WfIm1IVm = 2.0 * k1 * k2 * (1.0 - 3.0 * kappa) * kappa * pow((1.0 - 3.0 * kappa) * IVm + kappa * Im - 1.0, 2.0) * exp_1 + k1 * (1.0 - 3.0 * kappa) * kappa * exp_1;
        
        WfIVmIm1 = 2.0 * k1 * k2 * (1.0 - 3.0 * kappa) * kappa * pow((1.0 - 3.0 * kappa) * IVm + kappa * Im - 1.0, 2.0) * exp_1 + k1 * (1.0 - 3.0 * kappa) * kappa * exp_1;
        
        WfIVmIVm = 2.0 * k1 * k2 * pow(1.0 - 3.0 * kappa, 2.0) * pow((1.0 - 3.0 * kappa) * IVm + kappa * Im - 1.0, 2.0) * exp_1 + k1 * pow(1.0 - 3.0 * kappa, 2.0) * exp_1;
        
        // cout << "WfIm1Im1= " << WfIm1Im1 << endl;
    }
    
    if(VIm > 1.0)
    {
    	double exp_2 = exp( k2 * pow((1.0 - 3.0 * kappa) * VIm + kappa * Im - 1.0, 2.0));

        WfIm2 = k1 * kappa * (kappa * Im + (1.0 - 3.0 * kappa) * VIm - 1.0) * exp_2;
        WfVIm = k1 * (1.0 - 3.0 * kappa) * (kappa * Im + (1.0 - 3.0 * kappa) * VIm - 1.0) * exp_2;
        
        // Calculation of second derivatives
        
        WfIm2Im2 = 2.0 * k1 * k2 * pow(kappa,2.0) * pow((1.0 - 3.0 * kappa) * VIm + kappa * Im - 1.0, 2.0) * exp_2 + k1 * pow(kappa, 2.0) * exp_2;
        
        WfIm2VIm = 2.0 * k1 * k2 * (1.0 - 3.0 * kappa) * kappa * pow((1.0 - 3.0 * kappa) * VIm + kappa * Im - 1.0, 2.0) * exp_2 + k1 * (1.0 - 3.0 * kappa) * kappa * exp_2;
        
        WfVImIm2 = 2.0 * k1 * k2 * (1.0 - 3.0 * kappa) * kappa * pow((1.0 - 3.0 * kappa) * VIm + kappa * Im - 1.0, 2.0) * exp_2 + k1 * (1.0 - 3.0 * kappa) * kappa * exp_2;
        
        WfVImVIm = 2.0 * k1 * k2 * pow(1.0 - 3.0 * kappa, 2.0) * pow((1.0 - 3.0 * kappa) * VIm + kappa * Im - 1.0, 2.0) * exp_2 + k1 * pow(1.0 - 3.0 * kappa, 2.0) * exp_2;
    }
}

// ***

void Composite::calcW( const VectorXd & C )
{
	 matGround->calcW(C);
	 W = matGround->getW();

	 if (!hypertype.compare("coupled"))
     { 
         cout << "**ERROR* This model is to work with the decoupled form of the free-energy function" << endl;
     }
     else
     {
         calcAnisoTerms(C);
         
         if( IVm >= 1.0 and VIm >= 1.0 )
         {
              Wf = k1 / k2 / 2.0 * ( exp( k2 * pow( kappa * Im + (1.0 - 3.0 * kappa) * IVm - 1.0, 2.0)) - 1.0) + k1 / k2 / 2.0 * ( exp ( k2 * pow( kappa * Im + (1.0 - 3.0 * kappa) * VIm - 1.0, 2.0)) - 1.0);
         }
         else if ( IVm >= 1.0 and VIm < 1.0 )
         {
              Wf = k1 / k2 / 2.0 * ( exp( k2 * pow( kappa * Im + (1.0 - 3.0 * kappa) * IVm - 1.0, 2.0)) - 1.0);
         }
         else if (IVm < 1.0 and VIm >= 1.0 )
         {
              Wf = k1 / k2 / 2.0 * ( exp( k2 * pow( kappa * Im + (1.0 - 3.0 * kappa) * VIm - 1.0, 2.0)) - 1.0);
         }
         else
         {
        	  Wf = 0.0;
         }
     }
	
      W += Wf;
}

// ***

void Composite::calcS( const VectorXd & C )
{
    if (!hypertype.compare("decoupled"))
    {
    	matGround->calcS(C);	// Needs to be right here so that the ground invariants are properly computed.
        
        calcAnisoTerms(C);
        
        //  cout << "IVm (S)= " << IVm << endl;
        
        if (IVm >= 1.0 and VIm >= 1.0)
        {
            mSfTens = 2.0 * ((WfIm1 + WfIm2) * Id + WfIVm * A0 + WfVIm * B0);
        }
        else if (IVm >= 1.0 and VIm < 1.0 )
        {
            mSfTens = 2.0 * (WfIm1 * Id + WfIVm * A0);
            //cout << "Scomp= " << endl << mSfTens << endl;
        }
        else if (IVm < 1.0 and VIm >= 1.0 )
        {
            mSfTens = 2.0 * (WfIm2 * Id + WfVIm * B0);
        }
        else
        {
            mSfTens = 0.0;
        }
        
        matGround->calcP(C);
        SfTens = pow(matGround->getJ(), -2.0/3.0) * doubleContraction(matGround->getP(), mSfTens);
        
        VoigtFormat(SfTens, Sf);
    }
    else
    {
        cout << "**ERROR* This model is to work with the decoupled form of the free-energy function" << endl;
    }

    Svol = matGround->getSvol();
    Siso = matGround->getSiso() + Sf;
    S = matGround->getS() + Sf;

    //cout << "S= " << endl << S << endl;
}

// ***

void Composite::calcSiso(const VectorXd & C)
{
	calcS(C);
}

// ***

void Composite::calcCt( const VectorXd & C)
{
    if (!hypertype.compare("decoupled"))
    {
        matGround->calcCt(C); // Needs to be right here so that the invariants are properly computed.

        calcAnisoTerms(C);

        //cout << "IVm (Ct)= " << IVm << endl;
        //cout << "VIm (Ct)= " << VIm << endl;
        
        if (IVm >= 1.0 and VIm >= 1.0)
        {
            mCtfTens = 4.0 * pow(matGround->getJ(), -4.0/3.0) * ((WfIm1Im1 + WfIm2Im2) * selfDyadic(Id) + WfIVmIm1 * dyadic(Id,A0) + WfVImIm2 * dyadic(Id,B0) + WfIm1IVm * dyadic(A0,Id) + WfIm2VIm * dyadic(B0,Id) + WfIVmIVm * selfDyadic(A0)  + WfVImVIm * selfDyadic(B0));
        }
        else if (IVm >= 1.0 and VIm < 1.0)
        {
            mCtfTens = 4.0 * pow(matGround->getJ(), -4.0/3.0) * (WfIm1Im1 * selfDyadic(Id) + WfIVmIm1 * dyadic(Id,A0) + WfIm1IVm * dyadic(A0,Id) + WfIVmIVm * selfDyadic(A0));
        }
        else if (IVm < 1.0 and VIm >= 1.0)
        {
            mCtfTens = 4.0 * pow(matGround->getJ(), -4.0/3.0) * (WfIm2Im2 * selfDyadic(Id) + WfVImIm2 * dyadic(Id,B0) + WfIm2VIm * dyadic(B0,Id) + WfVImVIm * selfDyadic(B0));
        }
        else
        {
            mCtfTens = 0.0;
        }
    }
    else
    {
        cout << "**ERROR* This model is to work with the decoupled form of the free-energy function" << endl;
    }

    //matGround->calcCTensor(C);
    //matGround->calcCiTensor(C);
    matGround->calcP(C);
    matGround->calcPt(C);
    matGround->calctildeP(C);
    
    CtfTens = doubleContraction(doubleContraction(matGround->getP(), mCtfTens), matGround->getPt());
    CtfTens -= 2.0 / 3.0 * (dyadic(SfTens, matGround->getCiTensor()) + dyadic(matGround->getCiTensor(), SfTens));
    CtfTens += 2.0 / 3.0 * pow(matGround->getJ(), -2.0/3.0) * doubleContraction(matGround->getCTensor(),mSfTens) * matGround->gettildeP();

    VoigtFormat( CtfTens , Ctf);

   // cout << "CtTens= " << endl << CtTens << endl;
    
    Ctvol = matGround->getCtvol();
    Ctiso = matGround->getCtiso() + Ctf;

    Ct = Ctvol + Ctiso;
}

// ***

int Composite::getNfibers()
{	return nf;}

// ***
